In [61]:
import numpy
import matplotlib
import numpy as np
import matplotlib.cm as cm
import matplotlib.mlab as mlab
import matplotlib.pyplot as plt
%matplotlib inline
from pylab import *
Example: Two- dimensional Heat conduction -- Finite difference method & Lattice Boltzmann method
In this example we will use Adiabatic boundary condition for top and bottom boundary and hot and cold for left and right boundaries respectively. In this section we want to expand our code to two dimensional problems to become closed to module 4 which we learned in our numerical class.
What changes need to be considered?
1- For 2D we can use two different lattice positions D2Q5 and D2Q9. For heat conduction we can use D2Q5 however when we want to simulate Advection term D2Q9 necessarily gives us more accurate solution rather than D2Q5.
2- Notice to the defined lattice directions and their weighting coefficients.
3- No change will be made with collision unless we need to consider 4 other directions instead of just 2 (in 1D).
4- Please make attention how we create streaming procedure for two dimensional problem. If you make a mistake in streaming section then you will receive wrong answer (as like me).
In [62]:
nx=40 # the number of nodes in x direction lattice direction
ny=40 # the number of nodes in y direction lattice direction
alpha=0.25 # heat diffusion coefficient
D=2.0 # the dimension of the problem
mstep=2000 # the number of time step
omega=1./(0.5+(alpha*D)) #Chapman-Enskog dt=1 and dx=1
Tleft=1.0 # left wall temperature
Tright=0.0 # right wall temperature
k=5 # k=0,1,2,3,4,5,6,8,9
In [63]:
x=numpy.linspace(0,1,nx+1)
y=numpy.linspace(0,1,ny+1)
w=numpy.ones(k) # witghting factor
T=numpy.ones((ny+1,nx+1) ) # Temperature matrix
f= numpy.ones((k, ny+1,nx+1)) # distribution function
In [64]:
##================ Initial boundary condition
w[0]=0.0
w[1:5]=1./4.
In [65]:
##================== Initial value
T[0:ny+1,0:nx+1]=0.0
T[0:ny+1,0]=1.0
T[0:ny+1,nx]=0.0
T[ny,1:nx]=0.0
for i in range(nx+1):
for j in range(ny+1):
for l in range (k): #k=0,1,2,3,4
f[l,j,i]=w[l]*T[j,i]
In [66]:
## Main loop : comprised two parts :collision and streaming
##=====================
for n in range(mstep) :
# collision process
for i in range(nx+1):
for j in range(ny+1):
for l in range (k):
feq=w[l]*T[j,i]
f[l,j,i]=(1.-omega)*f[l,j,i]+omega*feq
#streaming process
# ==========================
for i in range(nx):
for j in range(ny,0,-1):
f[2,j,i]=f[2,j-1,i]
for i in range(nx,0,-1):
for j in range(ny,0,-1):
f[1,j,i]=f[1,j,i-1]
for i in range(nx,0,-1):
for j in range(0,ny):
f[4,j,i]=f[4,j+1,i]
for i in range(0,nx):
for j in range(0,ny):
f[3,j,i]=f[3,j,i+1]
# Boundary condition
# =============================
for j in range(0,ny+1) : #left Boundary
f[1,j,0]=( Tleft*(w[1]+w[3]) )-f[3,j,0]
for j in range(0,ny+1): #right Boundary
f[3,j,nx]=( Tright*(w[1]+w[3]) )-f[1,j,nx]
for i in range(0,nx+1): # bottom and top Boundary
f[4,ny,i]=-f[2,ny,i]
f[1,0,i]=f[1,1,i]
f[2,0,i]=f[2,1,i]
f[3,0,i]=f[3,1,i]
f[4,0,i]=f[4,1,i]
#================================
for i in range(nx+1):
for j in range(ny+1):
sum=0.0
for l in range (k):
sum=sum+f[l,j,i]
T[j,i]=sum
T[0:ny+1,0]=Tleft
T[0:ny+1,nx]=Tright
T[ny,1:nx]=0.0
T[0,1:nx]=T[1,1:nx]
#==============================
In [67]:
# Finite Difference Method
#=============================
Tf=numpy.ones((ny+1,nx+1) ) # finite difference Temperature matrix
To=numpy.ones((ny+1,nx+1) )
mstep=20000
dx=1./nx
dy=1./ny
alpha=0.25
dt=1e-4
Tf[0:nx+1,0:ny+1]=0.0
Tf[0:ny+1,0]=1.0
Tf[0:ny+1,nx]=0.0
Tf[ny,1:nx]=0.0
for n in range (mstep):
To[0:nx+1,0:ny+1]=Tf[0:nx+1,0:ny+1]
for j in range (1,ny):
for i in range (1,nx):
Txx=(To[j,i+1]+To[j,i-1])
Tyy=(To[j+1,i]+To[j-1,i])
Tf[j,i]=To[j,i]+ (dt*alpha*(Txx+Tyy-4.*To[j,i])/(dx**2))
#
Tf[0,1:nx]=Tf[1,1:nx]
##
In [68]:
CS1 = plt.contour(x,y,Tf)
title('Finite Difference Method')
clabel(CS1,inline='True')
plt.xlabel('$x$')
plt.ylabel('$y$')
show()
In [69]:
CS2 = plt.contour(x,y,T)
title('Lattice Boltzmann Method')
clabel(CS2,inline='True')
plt.xlabel('$x$')
plt.ylabel('$y$')
show()
In [70]:
plt.xlim(-0.001,1.1)
plt.ylim(-0.001,1.1)
T1=T[ny/2,0:nx+1]
T2=Tf[ny/2,0:nx+1]
plt.plot(x,T1, 'r-');
plt.plot(x,T2, 'bo');
plt.legend(['LBM','FDM']);
plt.xlabel('$x$')
plt.ylabel('$y$')
show()
In [71]:
mx, my = numpy.meshgrid(x,y)
plt.figure(figsize=(8,5))
plt.contourf(x,y,T,20)
plt.xlabel('$x$')
plt.ylabel('$y$')
plt.colorbar();